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ABSTRACT 

We introduce a numerical method to integrate tidal effects on collisional systems, using any 
definition of the external potential as a function of space and time. Rather than using a lin¬ 
earisation of the tidal field, this new method follows a differential technique to numerically 
evaluate the tidal acceleration and its time derivative. Theses are then used to integrate the mo¬ 
tions of the components of the collisional systems, like stars in star clusters, using a predictor- 
corrector scheme. The versatility of this approach allows the study of star clusters, including 
their tidal tails, in complex, multi-components , time -evolving external potentials. The method 
is implemented in the code NBODY6 dAarse thll2003l) . 
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1 INTRODUCTION 


Recent developments in algorithms and hardware open new per¬ 
spectives in the treatment of the collisional IV-body problem. 
For example, it is now possible to model the long-term evolu- 
tion of massive globular clusters with u p to several 10 5 stars 
dHeggiell201ll . l2014l : lHurlev & Sharall2012l) . Combined with much 
faste r methods allowing a wide exploration of the parameter space 


taster metnoas allowing a wide exploration ot me parameter space 
fe.g. Ijoshi et al.1 2000. Alexander & Gieles 2012, Hypki & Gierszl 


^ _ VDKl & Gl 

120131 Sollima & Mastrobuono Battistil 20141 Vasilievl 1 201 5| ), N- 
body simulations are reaching a high degree of realism, and 
can even be used to numerically reproduce real globular clusters 
dHeggiell2014l) . 

The old age of globular clusters usually requires to run sim¬ 
ulations covering an entire Hubble time to probe their full evolu¬ 
tion. However, cosmological simulations tell us that, during such 
a period of time, the galactic environments of these clusters un¬ 
dergo secular (accretion) and violent (galaxy interactions and merg¬ 
ers) evolutions (among many others, see e.g. lAgertz et alj|201 II) . 
The complexity of these external potentials is often neglected in 
TV-body simulations of clusters. However, several studies have 
adressed this problem, either by arbitrarily switching tidal effects 
to mimic the accretion of a dwarf satellite o nto a massive galaxy 
dMiholics et al.l 120141 : iBianchini et alJ I2015I) . or by (partly) cou¬ 
pling galaxy simulations to star cluster simul ations dFuiii et al.l 
120011: iRenaud & Gielesl2013l:lRieder et ai]|2013l) _ 

Among other approaches, in Rena ud et all d201 ll) . we pro- 
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posed a method to extract the tidal information (in the form of 
tables of tensors) along an orbit in a galaxy or cosmology sim¬ 
ulation. This method allows any kind of galactic potential (and 
cluster orbit) to be considered, including complex time-dependent 
ones like those found in galaxy mergers (see an application in 
iRenaud & Gieles! I20I3I) . It however suffers from three main limi¬ 
tations. 

(i) The need to run a galaxy simulation first. This can be numer¬ 
ically costly, time consuming and thus crippling for some users. 
Furthermore, running a full galaxy simulation is unnecessary when 
the evolution of the galaxy can be described analytically. (For ex¬ 
ample, the secular mass growth of a galaxy can be mimicked, at first 
order, by scaling the total mass and scale radii of the galaxy with- 
out modifying the s hape of its potential, see iDiemer et alJ[20131 : 
iBnist & Helmil2014l .I 

(ii) The different timescales between galactic and cluster scales. 
The tidal information is generally sampled at a much coarser fre¬ 
quency (~ 1-5 Myr and ~ 10-50 pc) than what is required in 
cluster simulations. Therefore the evaluation of this information at 
the time and position where it is needed by the cluster simulation 
requires interpolations of the data available. 

(iii) The so-called tidal approximation, i.e. the linearisation of 
the tidal forces. This introduces errors at large distances from the 
cluster centre and thus forbids the study of tidal tails. 

In this paper, we propose an alternative version which over¬ 
comes these limitations. In this new method, the user provides a 
program routine returning the galactic potential as a function of 
position and time. (Running a galaxy simulation beforehand is not 
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required, and the information on the external potential is not dis- 
cretised, which circumvents the need for interpolations.) The code 
uses this routine to compute the motion of the cluster, the tidal ac¬ 
celeration on its stars and the relevant derivatives (used to increase 
the accuracy, as described in Section 0. Since the galactic poten¬ 
tial is known at all possible positions, the tidal acceleration can be 
added to the motion of every star in the simulation, including those 
in tidal debris, which eliminates the limitation of the tidal approxi¬ 
mation. 

In the present version, the external potential must be defined 
by the user, in a form of a numerical code routine. It can be an sim¬ 
ple analytical function of position and time, and/or involve the nu¬ 
merical solving of more complex functional forms. This allows for 
a wide diversity of cases, from spherically symmetric static mod¬ 
els to time-evolving, triaxial, multi-component models. In some 
cases however, describing the potential with a function is too in¬ 
volved (e.g. galaxy me rgers), and it would be preferable to follow 
the approach of iRenaud et alj d201 ll) . i.e. tables of tensor coeffi¬ 
cients. Both methods have been implemented in NB0D Y6 and its 
GPU version, l lAarsetbir2003l ; lNitadori & Aarsethll 20121) . under the 
name NBODY 6tt, and are available on lineal 


2 METHOD 

In the rest of the paper, we describe the new method in the context 
of modelling a star cluster in a galactic potential. The method can 
however be used in other configurations, when an external potential 
must be included in a collisional iV-body system. 


2.1 Direct or differential? 

The evaluation of the contribution of the galaxy on a star cluster 
can be done in two ways. 

(i) Direct approach: the galactic acceleration on a star is added 
to that from the N — 1 other stars of the cluster. In that case, the 
coordinate frame is centered on the galaxy. 

(ii) Differential approach: the contribution of the galaxy is eval¬ 
uated at the position of the star and at the position the cluster (usu¬ 
ally its center of mass), and the difference is computed. The motion 
of the star is thus integrated with respect to the cluster, which has 
its own motion around the galaxy. The differential terms are called 
“tidal”. 

Because the galactic contribution on a star (e.g. on its accelera¬ 
tion) can be several orders of magnitude different than that from the 
cluster, the first method which consists in summing the two contri¬ 
butions can lead to numerical errors. For this reason, we adopt the 
second approach (as in all versions of NBODY 6 and NBODY 6tt). 


2.2 Numerical derivatives 

When the galactic potential yields an analytical description, the 
most accurate approach to include tidal effects would be to derive 
the galactic acceleration (and higher order terms) analytically. In 
some cases however, these derivations could be rather involved and 


1 http://personal.ph.surrey.ac.uk/~fr0005/nbody6tt.php 


even dissuasive. Such situations are encountered for complex func¬ 
tional forms of the potential (asymmetric, time-dependent etc), or 
in the rarer cases when the potential (or its derivatives) requires a 
numerical evaluation. If accurate enough, a numerical approach to 
compute both the orbit of the cluster and the tidal acceleration on its 
stars is often more appropriate and/or convenient. In the following 
Sections, we propose such a method and evaluate the error intro¬ 
duced. We found that, when solving the iV-body problem on single 
precision hardware (e.g. on Graphical Processing Unit, GPU) as it 
is often the case, the error introduced by our method gets truncated, 
such that our approach is as accurate as an analytical derivation, at 
least for the cases considered below (see Section[3}. 

In our method, a routine returns the galactic potential rp G 
as a function of position r and time t. As in the 2011 version 
of NBODY6tt, all computations are done in an inertial reference 
frame, where there are no fictitious forces (centrifugal, Coriolis and 
Euler, which are a priori unknown in the general case). The acceler¬ 
ation (per unit of mass) from the galaxy at the position r = ^ xtei 
and time t is thus given by minus the gradient of the potential, i.e. 


a G (r,t) 


d<p G (r,t) 

dr 


( 1 ) 


Numerically, we evaluate its i-th component using a central finite 
difference (see Fig. |TJ with a fourth-order accuracy: 


a G (r,t) tts — (12 hi) 1 [(pair — 2hiei,t ) — 8rp G (r — hid,t) 
+8<pa(r + hi€i, t ) — cp G (r + 2 hid, f)] , (2) 

which thus has errors O(hf). The choice of the value of the step 
size hi is described in Section [2~5l 

The time derivative of the galactic acceleration, or jerk (jo), 
is then computed using the relation 


ja{r,t) = 


da a (r,t) 

dt 


da G (r,t) dr da G (r,t) 


dr 

= T (r,t)v + 


dt dt 
da G (r,t) 
dt ’ 


(3) 


where v is the velocity vector of the star with respect to the galaxy. 
The space derivative of the galactic acceleration (i.e. minus the 
Hessian matrix of the galactic potential) is the tidal tensor T, which 
we compute numerically by using a second order central finite dif¬ 
ference of <p G , with a second-order accuracy (Fig. 0 . Its compo¬ 
nents read 


T xj (r,t) ft* —(4 hihj)- 1 {(pair — hid + hjej,t) 

—<p G {r + hid + hjej,t ) 

+(pa(r + hiei - hjej,t) 

-<po(r - hid - hjej,t) ]. (4) 

The partial time derivative of the acceleration, computed with a first 
order finite difference of a G with a second-order accuracy, is 

daa ^ ,t ) ~(2/r i ) —1 [—a G (r, t-h t ) + a G {r, t + h t )] (5) 

where ht is the step size for the time dimension. Here, the fourth- 
order accuracy provided by equation 0 is unnecessary. Instead, we 
compute the acceleration with a second-order accuracy: 

a G (r, t) « - (2hi)~ 1 [—cp G (r - hid,t) + cp G (r + hid,t)] , 

(6) 
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Figure 1. Stencil used to compute the first (acceleration) and second (jerk) 
order derivatives of the galactic potential at a given position (green star). 

such that the z-th component of its partial time derivative reads 

d a °ir ^ ~ [ < fa( r ~ hi e i, t - ht) 

~4>g {r + hiei,t — h t ) 

—(pair — hiei,t + ht) 

+(p G (r + hiei,t +h t )}. (7) 

For simplicity, we use the same values of the spacial step sizes {hi} 
at t, t — ht and t + ht- Finally, the jerk is evaluated by replacing 
equations 0 and 0 in equation 0. 

2.3 Cluster orbit 

The motion of the cluster around the galaxy is described using a 
guiding centre, as already done in NBODY6 by setting the option 
KZ(14) to 3 or 4. This pseudo-particle initially matches the cen¬ 
tre of mass of the cluster, but can slightly deviate from it later 
on, as stars tidally ejected from the cluster take away momen¬ 
tum in an asymmetric fashion. NB0DY6tt integrates the equa¬ 
tion of motion of the guiding center using the galactic acceleration 
(equation 0 and jerk (equation 0 in a predictor-corrector scheme 
lAarsethll20031) . 

Note that dynamical friction is not included in the integration 
of the cluster orbit (but see Petts et al., in preparation). 


2.4 Tidal acceleration 

When the regular force on a star (either in the cluster or in the tidal 
debris) must be updated, the contribution of the galaxy is evalu¬ 
ated, using the differential approach presented in Section [2711 The 
galactic acceleration (resp. jerk) at the position of the guiding cen¬ 
ter of the cluster is subtracted from that at the position of the star. 
We thus obtain the tidal (i.e. differential) acceleration (resp. jerk), 
used in the predictor-corrector integration scheme of NB0DY6 to 
evolve the motion of the star in the cluster. The advantage of this 
approach is that it allows us to follow the motion of the stars on 
galactic orbits after they have left the cluster. 


2.5 Step sizes 

In this Section, for simplicity, we consider the derivation of the 
potential rpo with respect to a single dimension x, with the step size 
h, and we will generalize our approach to the multi-dimensional 
case later. 

The accuracy in the evaluations of the numerical derivatives 
presented in Section02]rely on the choice of step sizes. A too small 
value would lead to round-off errors while the derivative would not 
be accurate for a too large value. According to lPress et al.1 d2007L 
their Section 5.7), the optimum value of h depends on the “curva¬ 
ture scale” x c of the potential at the position it is evaluated, as 

h ~ e^Xc, (8) 

where e is comparable to the machine accuracy (i.e. ~ 10” 16 
in double precision). Both £ and x c depend on the order of the 
derivative considered. By writing the Taylor series expansion of 
4>g{x + h) and by seeking the value of h which minimizes the sum 
of the round-off and truncation errors (see an example at lower or¬ 
der in lPress et al.tt2007t their Equation 5.7.5), one can show that 

/ \ 1/5 

and < 9 > 

\ dx b / 

for equation 0. and 

/ \ 1/3 

< = l and = bfe ( |0 > 

\ dx 3 / 

for equations 0 and 0. Therefore, computing the optimum step 
size requires, paradoxically, the evaluation of higher order deriva¬ 
tives of the potential. In principle, this could be done through the 
approach of lRiddersI d 19821) . at the cost of many additional evalua¬ 
tions of the potential. Once the optimum step size found, we could 
build the stencil of Fig0and compute all the relevant derivatives, 
for all stars in the cluster. In practice, this methods appears to be of 
limited interest with resp ect to its significa nt extra numerical cost. 
Following the advice of IPress et alj d2007h . we choose to adopt a 
much simpler and faster approach and assume that the curvature 
scale can be approximated by the position x. In that case, Fig. 0 
shows the relative numerical error made on the acceleration and 
the jerk using the potential of a point-mass (tpa = — 1 / at). To take 
avantage of the stencil method (Fig.0, and to further limit the num¬ 
ber of evaluations of the potential, we choose to use the same step 
size for both the computation of the acceleration and the jerk. Since 
the jerk is only used in the predictor corrector method, the accuracy 


© 0000 RAS, MNRAS 000, 000-000 














































4 Renaud & Gieles 



Figure 2. Relative numerical error f 1 - numerical/exact) made on the accel¬ 
eration (red) and the jerk (blue) from a point-mass potential, as a function 
of the step size (for x = 1). The vertical dotted line marks the choice of h 
adopted in our implementation (equation ! 1 It . 


an “hardware truncation” of the accuracy of the numerical variables 
they manipulate. Such truncation affects the evaluation of the reg¬ 
ular force in NB0DY6, to which the tidal force is added. Therefore, 
despite being evaluated in double-precision, the tidal force applied 
to the particles are truncated to single-precision when using GPUs. 
This limitation does not concern simulations run on CPUs. 


2.7 Energy conservation 

One way of monitoring the numerical errors made during a sim¬ 
ulation is to control the conservation of energy. However, time- 
dependent tides imply that the energy is not conserved. The code 
circumvented this issue through the following methotJE Let W be 
the time derivative of the internal energy of the cluster (potential 
energy U from internal interactions plus kinetic energy K relative 
to the guiding centre of the cluster). It can be written as 


w= d{u + K) 

dt 



d U dr; 
dr; df 


+E 


rriiVi- 


dvj 

dt 


N N 

= ^2 —miaciVi + ^2 miV^aci + aa) 
i i 


on its value is less critical than that of the acceleratioiQ. Therefore, 
we have adopted the optimum step size for the acceleration to com¬ 
pute both the acceleration and the jerk. 

Setting h oc x implies numerical issues for x = 0 leading 
to infinite derivatives. To limit the number of occurrences of this 
situation, we choose instead to use h oc r. The step size would then 
be far from being optimum where, e.g. x -C r, but the error would 
generally be made on a small component of the total acceleration 
(a G 11 qq 11). In the end, we adopt: 

ft; = 4 x 10' 4 r (11) 

for all three values of i, and use this empirical relation for all po¬ 
tentials. We note however that our choice of the step size would 
not be optimum at the vicinity of substructures in the potential, like 
spirals arm in a galactic disc. 

For time-dependent potentials, the time step ht is taken to be 
the same as the time step of the predictor corrector scheme used 
when integrating the motion of the guiding centre in the galaxy. It 
is the same for all stars. 

This method is tested in Section[3] 


2.6 Numerical precision 

In the implementation of our method in NBODY 6, the computations 
of the external potential and its derivatives are performed in double¬ 
precision to minimize the impact of the loss of accuracy during the 
numerical derivations. 

We note however that the actual force might be truncated to 
single-precision in some cases. The majority of the GPUs used in 
the community are limited to single-precision and thus introduce 


N 

= UliViCLGi, (12) 

i 

where ac represents the internal acceleration due to the stars in the 
cluster, rrii is the mass of the i-th star, and the sums are made over 
all stars in the system. We note that the second time derivative of 
W reads 

N 

W = ^2 rrii [(ac; + a G ;)a G ; + Vij Gi ] • (13) 


The definition of W implies that 

U + K — J W dt = constant, 


(14) 


which is the assertion that the code must verify. 

In practice, to numerically obtain the variation ( W) of internal 
energy between the timesteps ti and to, we compute the Taylor 
series of W at the time to, truncated to the second order: 


W(ti) 


dW 

w (t 0 ) + ^ r (t 1 -to) + 


d 2 W (ti — to) 2 
dt 2 2 


(15) 


whence 


AW &WAt + W^-. (16) 

Using equations 01 and O, we can compute the variation AW 
at a given timestep in the simulation. The accumulation of these 
variations since the beginning of the simulation gives the numerical 
equivalent of f W dt that we subtract to the value of U + K to 
check equation ( Il4pl 


2 The jerk is multiplied by the small time step of the predictor-corrector 
scheme before being added to the acceleration. Therefore, it is the error on 
the jerk times this time step (typically 10 —6 the orbital period) which must 
be compared to the error on the acceleration, such that the dominant source 
of errors is that made on the acceleration. 


3 already used in previous versions of NBODY 6 when setting the option 
KZ (14) =3. 

4 From equation ®. one can see that — f W dt = — W is the part of the 
tidal energy that contributes to energy conservation. It is not the full tidal 
energy in the general case. 
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Figure 3. Relative numerical error (1 - numerical/exact) on the x- 
component of the acceleration, as a function of x. The potentials adopted 
have the form —1/x (black) and —(1 + .r 2 ) ~ 1 / 2 (red). The error remains 
small over the range considered, despite an increase when the potential be¬ 
comes flatter and flatter. 

3 TESTS 

In this Section, we compare the results from our implementation of 
the method to either analytical solutions or numerical results from 
NB0DY6, and thus do not explore the full possibilities offered by 
the new method. 


3.1 Acceleration 

Fig-HJshows the relative numerical error made on the ^-component 
of the acceleration (evaluated with equation^. We monitor the ac¬ 
curacy of the derivation scheme by both setting divergent (cusped) 
and cored potentials, such that a wide variety of potential slopes 
are considered. The relative error on the acceleration computed 
in double-precision overcomes 10“', i.e. the machine accuracy in 
single-precision, when it is computed at a distance to the centre 
of the potential ~ 3 x 10 6 times shorter than the potential char¬ 
acteristic scale. In the context of cored potentials of galaxies, this 
would corresponds to cluster-galaxy distances of a few parsecs, i.e. 
a relatively rare situation. 

Therefore, except in the case previously mentioned, following 
the discussion in Section [231 the net accuracy of the acceleration 
would be set by the machine precision when using single-precision 
GPUs. In that view, the acceleration evaluated with our method 
could be considered as accurate as a computation using the ana¬ 
lytical expression. 

3.2 Cluster orbit 

Fig. |4] shows the numerical errors made in the energy and angular 
momentum of the guiding center of a cluster on an eccentric or¬ 
bit of eccentricity 0.5 and apocenter distance 3 kpc around a point 
mass galaxy (10 9 M@). This corresponds to an orbital period of 
~ 260 Myr. The errors from the new method remain of the or¬ 
der of those from the original method. A long-term drift exists but 



Figure 4. Relative numerical errors (1 - numerical/exact) on the orbital 
energy (top) and angular momentum (bottom), of the guiding centre of a 
cluster on an eccentric orbit (eccentricity 0.5, apocenter = 3 kpc) around a 
point-mass galaxy (10® Mg). 


is limited to about one dex over 10 orbital periods, indicating that 
both quantities are well conserved over periods of time matching 
the typical life-time of clusters, or the typical duration of these sim¬ 
ulations. 


3.3 Cluster evolution 

Fig.[5]compares the evolution of the mass and a few Lagrange radii 
of clusters modelled with the new method with that of the same 
clusters modelled with NBODY6. The clusters are equal-mass mod¬ 
els of 4096 or 8192 stars (1 Mq each) on a Plummer profile with 
an initial virial radius of 1 pc, and placed on an eccentric orbit 
around a point-mass galaxy, as described in the previous Section. 
The quantities plotted are computed using the bound stars, i.e. those 
for which the sum of the kinetic and internal energy (with respect 
to the guiding centre) is negative, as in iRenaud et al.1 d20 1 ill . (In 
other words, we neglect the tidal energy when deciding the mem¬ 
bership of stars.) To evaluate the amplitude of Poisson’s noise in 
our measurements, we have run simulations of four realisations of 
each cluster by changing the seed of the random number generator 
used to produce the initial conditions. 

The evolution of the mass and size of the clusters is com¬ 
pared to that computed using the built-in version of NBODY6 (op¬ 
tion KZ ( 14 ) =3). We note that the low-frequency evolution of the 
clusters, mainly connected to the orbital period, is well reproduced 
by the new implementation until the end of the simulations. As ex¬ 
pected, the largest deviations between the two methods are found in 
the inner regions of the cluster, where the evolution is the most vul¬ 
nerable to statistical effects leading to the stochastic formation and 
destruction of binaries. The amplitude of the deviations increases 
as the number of stars shrinks, as a result of Poisson noise. 

In conclusion, the agreement between the two methods is very 
good for the entire lifetime of the clusters considered, and over all 
phases of their orbit (apocenter, pericenter and in between). 


© 0000 RAS. MNRAS 000, 000-000 
























6 Renaud & Gieles 



0 1000 2000 3000 4000 

time [Myr] 




0 2 4 6 8 10 12 14 

time [Myr] 


Figure 6. Top: evolution of the internal energy (potential + ki- 
netic), normalised t o its initial value, of a cluster plunging through a 
iMiyamoto & Nagail ll975t) disc of scale-heigh b = 300 pc (blue) and of 
b = 10 pc (red). The instants when the cluster reaches the scale-heigh of 
the disc (z = d=6) and the median plan (z = 0) are marked with crosses 
and a plus-sign. Bottom: relative difference in the internal energy between 
the two methods. The dashed curves correspond to the setup where the ini¬ 
tial velocity is given to the disc, and not the cluster (Section [3.51 . (Spikes 
in the curves corresponds to high velocity stars being ejected from the clus¬ 
ter. Such events are subject to Poison noise and thus varies from method to 
method.) 


Figure 5. Evolution of the mass and the structure (shown with the 10%, 
50% and 90% Lagrange radii) of models of a clusters with initially 4096 and 
8192 stars (labelled 4k and 8k respectively) integrated with the new method, 
and compared to the NBODY6 original method. Each curve represents the 
median value over four 7V-body realisations of the initial conditions, and 
the shaded areas show the standard deviations of these realisations (origi¬ 
nal method only), following lEmst et~akl 1201 ll) and Whitehead et al. |201 3| ). 
Strong variations of the radii near the time of dissolution of the clusters are 
due to uncertainties in the determination of the centre of small-iV systems. 


3.4 Steep potentials and tidal shocks 

To further test the method, we consider the case of a cluster plung¬ 
ing through a disc. The ext ernal potential is modelled with a single 
iMivamoto & Nagai] J 1975h disc: 


x 2 + y 2 + (a + d z 2 + b 2 ) 

with the parameters M = 10 11 Mg, a = 5 kpc and b = 300 pc. 
The cluster is initially set at the position xo = yo = 0, zq = 
1.5 kpc, i.e. above the median plane of the disc, with the ini¬ 
tial velocity toward the disc v z ,o = —150 km s -1 . Such config¬ 
uration can be setup in the original version of NBODY6 (option 
KZ (14 ) =3) and thus we can compare here again the result from 



the new method to that of NBODY6. The cluster is made of 8196 
equal-mass stars distributed on a Plummer profile with an initial 
virial radius of 3 pc. The cluster crosses the median plan of the 
disc with a velocity of « 230 kms' 1 and takes « 2.6 Myr to 
cover the 2 x 300 pc of the disc scale-height. 

Fig. [6] shows the tidal heating of the cluster, and the relative 
difference in the internal energy of the cluster (potential + kinetic) 
between the two methods. The encounter leads to an increase of 
about one percent of the cluster initial internal energy over a few 
Myr. During this period, the rapidly varying external potential in¬ 
duces differences results from two methods. It is likely that our 
choice of assuming a universal curvature scale (see Section [2~5t is 
responsible for most of the differences. However, the relative dif¬ 
ference remains small (~ 10 -4 ). 

To push the code to its limits, we consider the extreme case of 
a cluster moving through a very thin disc. We used the same poten¬ 
tial form as before, but with a scale-height b = 10 pc (i.e. about 
3 times the initial virial radius of the cluster). We set the cluster 
at zo = 1 kpc with the same initial velocity as before (v z ,o = 
— 150 km s _1 ). The “impact” velocity is ss 225 km s _1 , meaning 
that the cluster takes « 0.09 Myr to cover the 2 x 10 pc of the disc 
scale-heig ht, i.e. shorter t han t he cluster crossing t i me (~ 1 Myr). 
Following ISnitzerl dl987t) and iGnedin & OstrikeH d 19971) . we can 
consider such encounter as an impulsive tidal shock, as opposed to 
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an adiabatic effect. The comparison between the original and new 
methods is showed in Fig. [6] Despite an energy gain about twice 
larger than in the thicker disc case, the relative difference between 
the two method remains of the same amplitude as before. 

To conclude, the difference in internal energy of a cluster com¬ 
puted with the two methods remains well below one percent, even 
in steep potentials (i.e. with a curvature scale being a strong func¬ 
tion of position). 

3.5 Time-dependent potential 

Finally, we test the method in the context of a time-dependent po¬ 
tential. Note that, in our approach, time-dependence only affects 
the predictor-corrector scheme (thought the time derivative in the 
expression of the jerk, equation [7]l and not the (tidal) acceleration 
itself which is computed at each timestep in a static way, using the 
present-time expression of the potential. 

NBODY 6 does not allows for a treatment of tides with an ex¬ 
plicit time dependence. Therefore, to compare our method to the 
original code, we adopt the following approach. 

In the original NBODY 6 , we use the same setup as in the pre¬ 
vious Section, i.e. a cluster plunging through a static disc with an 
initial velocity v z ,o- In the new version (NBODY6tt), we setup the 
cluster with no initial velocity, and define the external potential as 
the same disc but moving toward the cluster with the velocity —v Z: o 
by replacing a with 2 — v z $t in equation G3. such that our poten¬ 
tial becomes time-dependent. The two setups are equivalent and the 
physical evolution of the system should be exactly the same in both 
cases. The relative difference in internal energy (for the two val¬ 
ues of b adopted before, i.e. 300 and 10 pc) is showed in Fig. [6] 
The differences induced by the disc-cluster encounter noted in the 
previous Section are still present here. Furthermore, the offset of 
the potential from the origin makes our choice of the step-size as a 
linear function of position (equation ! 1 11 1 less optimum than before. 
As a consequence, the differences between the two methods are 
slightly amplified than for the case of the static potential centered 
on the origin. Despite a larger amplitude than the static cases ex¬ 
plored before, the relative differences remains of the order of 1CP 4 , 
showing that our method has a comparable behaviour than the orig¬ 
inal NBODY 6 , even in such extreme cases. 


4 CONCLUSION 

We introduce a method to compute the tidal acceleration on an col- 
lisional JV-body system embedded in an external potential which 
can be described with a function of position and time (analytical 
and/or numerical). The method evaluates the first and second space 
derivatives of the potential to obtain the tidal acceleration and the 
tidal tensor. The tensor allows us to estimate the tidal jerk which, 
with the acceleration, is used in a predictor-corrector scheme to in¬ 
tegrate the equations of motion of the A-bodies. By circumventing 
the need of the classical tidal approximation (linearisation of the 
tidal force), this method can accurately integrate the motion of any 
body in the system, including those in tidal debris. 

The orbit of the guiding center of the system within the ex¬ 
ternal potential is also computed, following a comparable method. 
The numerical errors made on these quantities are of the order of 
10 -12 or smaller (for the test cases we considered), and thus the 


resulting net evolution of the A-body system is comparable to that 
obtained with other approaches. We have considered several setups 
where the evolution of clusters could be compared to that provided 
by existing methods, and found reasonable agreements in all cases. 
This suggests that the new method is able to produce simulations 
with accuracy standards close to that of NBODY 6. 

This new method however allows a larger flexibility, as any 
external potential can be considered, providing it can be described 
by a numerical routine. Among the endless list of possible applica¬ 
tions, one can imagine to describe the tidal effects of time-evolving 
multi-component galaxies including halo, bulge, disc(s), spiral pat- 
tern(s), bar(s), ring(s), and/or undergoing accretion of intergalactic 
gas and (to some extent) satellite galaxies. Alternative methods are 
however required when the external potential cannot be described 
by a numerical function, like for example in major galaxy mergers. 
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